collegeSummaryRTD

Предварительный анализ данных

library(readxl)
col_I_sn<-read_excel("stat2021_sm\\us_col.std\\I_shortname.xls")
## New names:
## • `` -> `...1`
col_I<-read_excel("stat2021_sm\\us_col.std\\I.xls")
## New names:
## • `` -> `...1`
head(col_I_sn)|>rmarkdown::paged_table()
summary(col_I_sn)
##      ...1               PPIND            FICE          STATE          
##  Length:176         Min.   :1.000   Min.   : 1009   Length:176        
##  Class :character   1st Qu.:1.000   1st Qu.: 1738   Class :character  
##  Mode  :character   Median :1.000   Median : 2550   Mode  :character  
##                     Mean   :1.301   Mean   : 2901                     
##                     3rd Qu.:2.000   3rd Qu.: 3417                     
##                     Max.   :2.000   Max.   :10366                     
##                                                                       
##      TYPE              AVRMATH         AVRVERB         AVRCOMB    
##  Length:176         Min.   :390.0   Min.   :391.0   Min.   : 810  
##  Class :character   1st Qu.:520.0   1st Qu.:454.5   1st Qu.: 980  
##  Mode  :character   Median :544.0   Median :479.5   Median :1017  
##                     Mean   :563.5   Mean   :494.7   Mean   :1058  
##                     3rd Qu.:603.0   3rd Qu.:518.5   3rd Qu.:1126  
##                     Max.   :750.0   Max.   :665.0   Max.   :1410  
##                     NA's   :68      NA's   :68      NA's   :67    
##     AVR_ACT          MATH_1          MATH_3          VERB_1     
##  Min.   :19.00   Min.   :350.0   Min.   :480.0   Min.   :320.0  
##  1st Qu.:22.00   1st Qu.:460.0   1st Qu.:590.0   1st Qu.:400.0  
##  Median :23.00   Median :500.0   Median :627.5   Median :430.0  
##  Mean   :23.69   Mean   :515.4   Mean   :635.0   Mean   :448.3  
##  3rd Qu.:25.00   3rd Qu.:570.0   3rd Qu.:687.5   3rd Qu.:480.0  
##  Max.   :31.00   Max.   :740.0   Max.   :780.0   Max.   :630.0  
##  NA's   :83      NA's   :34      NA's   :34      NA's   :34     
##      VERB_3          ACT_1           ACT_3          APP_REC     
##  Min.   :440.0   Min.   :16.00   Min.   :21.00   Min.   :  787  
##  1st Qu.:520.0   1st Qu.:19.00   1st Qu.:25.00   1st Qu.: 4323  
##  Median :550.0   Median :21.00   Median :26.00   Median : 7654  
##  Mean   :561.3   Mean   :21.28   Mean   :26.55   Mean   : 8544  
##  3rd Qu.:600.0   3rd Qu.:23.00   3rd Qu.:28.00   3rd Qu.:11776  
##  Max.   :720.0   Max.   :29.00   Max.   :33.00   Max.   :48094  
##  NA's   :34      NA's   :67      NA's   :67      NA's   :1      
##     APP_ACC         NEW_STUD        NEW10           NEW25          FULLTIME    
##  Min.   :  507   Min.   : 210   Min.   : 8.00   Min.   :24.00   Min.   :  912  
##  1st Qu.: 3033   1st Qu.:1264   1st Qu.:24.00   1st Qu.:52.00   1st Qu.: 5846  
##  Median : 4761   Median :1949   Median :32.00   Median :63.00   Median :10215  
##  Mean   : 5546   Mean   :2252   Mean   :41.48   Mean   :65.55   Mean   :11296  
##  3rd Qu.: 7232   3rd Qu.:3035   3rd Qu.:57.00   3rd Qu.:82.00   3rd Qu.:15033  
##  Max.   :26330   Max.   :7425   Max.   :98.00   Max.   :99.00   Max.   :31643  
##  NA's   :1                      NA's   :16      NA's   :26                     
##     PARTTIME          IN_STATE        OUT_STAT        R_B_COST   
##  Min.   :   16.0   Min.   :  647   Min.   : 2279   Min.   :2082  
##  1st Qu.:  804.8   1st Qu.: 2100   1st Qu.: 6712   1st Qu.:3588  
##  Median : 1694.5   Median : 3030   Median : 8400   Median :4213  
##  Mean   : 2519.9   Mean   : 6641   Mean   :10108   Mean   :4538  
##  3rd Qu.: 3240.8   3rd Qu.:12348   3rd Qu.:12668   3rd Qu.:5564  
##  Max.   :21836.0   Max.   :20100   Max.   :20100   Max.   :7425  
##  NA's   :6         NA's   :7       NA's   :1                     
##       ROOM          BOARD         ADD_FEE            BOOK           PERSONAL   
##  Min.   :1033   Min.   :1000   Min.   :  20.0   Min.   : 300.0   Min.   : 300  
##  1st Qu.:1810   1st Qu.:1762   1st Qu.: 210.0   1st Qu.: 500.0   1st Qu.:1164  
##  Median :2644   Median :2125   Median : 425.5   Median : 600.0   Median :1600  
##  Mean   :2708   Mean   :2241   Mean   : 648.1   Mean   : 603.1   Mean   :1763  
##  3rd Qu.:3490   3rd Qu.:2586   3rd Qu.: 694.0   3rd Qu.: 673.8   3rd Qu.:2200  
##  Max.   :6965   Max.   :4760   Max.   :4374.0   Max.   :1230.0   Max.   :6800  
##  NA's   :42     NA's   :64     NA's   :40       NA's   :2        NA's   :11    
##       PH_D           TERM_D         SF_RATIO         DONATE     
##  Min.   :63.00   Min.   :67.00   Min.   : 2.90   Min.   : 4.00  
##  1st Qu.:80.50   1st Qu.:87.00   1st Qu.:10.88   1st Qu.:10.75  
##  Median :87.00   Median :92.00   Median :14.50   Median :17.00  
##  Mean   :85.72   Mean   :90.52   Mean   :14.23   Mean   :19.01  
##  3rd Qu.:92.00   3rd Qu.:96.00   3rd Qu.:18.02   3rd Qu.:24.00  
##  Max.   :99.00   Max.   :99.00   Max.   :24.70   Max.   :54.00  
##  NA's   :9       NA's   :16                      NA's   :12     
##     INSTRUCT        GRADUAT         SAL_FULL          SAL_AC     
##  Min.   : 3605   Min.   :10.00   Min.   : 446.0   Min.   :364.0  
##  1st Qu.: 7604   1st Qu.:47.50   1st Qu.: 597.0   1st Qu.:445.0  
##  Median : 9840   Median :62.00   Median : 661.0   Median :479.5  
##  Mean   :12832   Mean   :62.02   Mean   : 669.2   Mean   :487.1  
##  3rd Qu.:14340   3rd Qu.:74.50   3rd Qu.: 732.2   3rd Qu.:521.2  
##  Max.   :62469   Max.   :99.00   Max.   :1009.0   Max.   :733.0  
##                  NA's   :5                                       
##      SAL_AS         SAL_ALL         COMP_FUL         COMP_AC     
##  Min.   :323.0   Min.   :362.0   Min.   : 537.0   Min.   :438.0  
##  1st Qu.:381.8   1st Qu.:472.2   1st Qu.: 729.0   1st Qu.:556.5  
##  Median :407.0   Median :522.5   Median : 815.0   Median :606.0  
##  Mean   :412.8   Mean   :534.0   Mean   : 827.2   Mean   :612.1  
##  3rd Qu.:434.2   3rd Qu.:578.2   3rd Qu.: 910.0   3rd Qu.:662.2  
##  Max.   :576.0   Max.   :866.0   Max.   :1236.0   Max.   :909.0  
##                                                                  
##     COMP_AS         COMP_ALL         NUM_FULL         NUM_AC     
##  Min.   :395.0   Min.   : 436.0   Min.   : 39.0   Min.   : 32.0  
##  1st Qu.:481.0   1st Qu.: 587.0   1st Qu.:184.5   1st Qu.:138.5  
##  Median :508.5   Median : 652.0   Median :278.0   Median :208.0  
##  Mean   :519.0   Mean   : 665.8   Mean   :336.8   Mean   :230.3  
##  3rd Qu.:552.2   3rd Qu.: 730.2   3rd Qu.:457.8   3rd Qu.:299.0  
##  Max.   :717.0   Max.   :1075.0   Max.   :997.0   Max.   :721.0  
##                                                   NA's   :1      
##      NUM_AS         NUM_INS          NUM_ALL      
##  Min.   : 29.0   Min.   :  0.00   Min.   : 108.0  
##  1st Qu.:128.5   1st Qu.:  5.00   1st Qu.: 505.5  
##  Median :175.0   Median : 16.00   Median : 721.0  
##  Mean   :190.7   Mean   : 27.09   Mean   : 812.8  
##  3rd Qu.:238.2   3rd Qu.: 40.00   3rd Qu.:1035.0  
##  Max.   :510.0   Max.   :178.00   Max.   :2261.0  
##                  NA's   :1
col_I_sn$PPIND=ifelse(col_I_sn$PPIND==1, "public", "private")

NA review

col_I_sn_noNA<-col_I_sn|>na.omit()
paste("col_I_sn dim: ", paste(dim(col_I_sn),collapse=", "))
## [1] "col_I_sn dim:  176, 49"
paste("col_I_sn_noNA dim: ", paste(dim(col_I_sn_noNA), collapse=", "))
## [1] "col_I_sn_noNA dim:  22, 49"
naCol<-colSums(is.na(col_I_sn))|>sort(decreasing=TRUE)
naCol
##  AVR_ACT  AVRMATH  AVRVERB  AVRCOMB    ACT_1    ACT_3    BOARD     ROOM 
##       83       68       68       67       67       67       64       42 
##  ADD_FEE   MATH_1   MATH_3   VERB_1   VERB_3    NEW25    NEW10   TERM_D 
##       40       34       34       34       34       26       16       16 
##   DONATE PERSONAL     PH_D IN_STATE PARTTIME  GRADUAT     BOOK  APP_REC 
##       12       11        9        7        6        5        2        1 
##  APP_ACC OUT_STAT   NUM_AC  NUM_INS     ...1    PPIND     FICE    STATE 
##        1        1        1        1        0        0        0        0 
##     TYPE NEW_STUD FULLTIME R_B_COST SF_RATIO INSTRUCT SAL_FULL   SAL_AC 
##        0        0        0        0        0        0        0        0 
##   SAL_AS  SAL_ALL COMP_FUL  COMP_AC  COMP_AS COMP_ALL NUM_FULL   NUM_AS 
##        0        0        0        0        0        0        0        0 
##  NUM_ALL 
##        0
table(naCol)
## naCol
##  0  1  2  5  6  7  9 11 12 16 26 34 40 42 64 67 68 83 
## 21  5  1  1  1  1  1  1  1  2  1  4  1  1  1  3  2  1

log dataset

q_cols<-sapply(col_I_sn, is.numeric)
q_cols[c("PPIND", "FICE")]<-FALSE
q_cols<-names(q_cols)[q_cols==TRUE]

col_I_sn_log<-col_I_sn|>
  mutate(across(q_cols, ~log(.x)))|>
  rename_with(~paste0("log-", .x), q_cols)

head(col_I_sn_log)|>rmarkdown::paged_table()

Разобраться в том, что означают признаки.

Выделим основные признаки, по которым будет проходить дальнейший анализ данных

Список признаков.

Количественные признаки:

  • AVRMATH Average Math SAT score
  • AVRVERB Average Verbal SAT score
  • AVRCOMB Average Combined SAT score
  • AVR_ACT Average ACT score
  • MATH_1 First quartile - Math SAT
  • MATH_3 Third quartile - Math SAT
  • VERB_1 First quartile - Verbal SAT
  • VERB_3 Third quartile - Verbal SAT
  • ACT_1 First quartile - ACT
  • ACT_3 Third quartile - ACT
  • APP_REC Number of applications received
  • APP_ACC Number of applicants accepted
  • NEW_STUD Number of new students enrolled
  • FULLTIME Number of fulltime undergraduates
  • PARTTIME Number of parttime undergraduates
  • IN_STATE In-state tuition
  • OUT_STAT Out-of-state tuition
  • R_B_COST Room and board costs
  • ROOM Room costs
  • BOARD Board costs
  • ADD_FEE Additional fees
  • BOOK Estimated book costs
  • PERSONAL Estimated personal spending
  • PH_D Pct. of faculty with Ph.D.’s
  • TERM_D Pct. of faculty with terminal degree
  • SAL_FULL Average salary - full professor
  • SAL_AC Average salary - associate professor
  • SAL_AS Average salary - assistant professor
  • SAL_ALL Average salary - all ranks
  • COMP_FUL Average compensation - full professor
  • COMP_AC Average compensation - associate professor
  • COMP_AS Average compensation - assistant professor
  • COMP_ALL Average compensation - all ranks
  • NUM_FULL Number of full professor
  • NUM_AC Number of associate professor
  • NUM_AS Number of assistant professor
  • NUM_INS Number of instructors
  • NUM_ALL Number of faculty - all ranks
  • INSTRUCT Instructional expenditure per student
  • GRADUAT Graduation rate
  • SF_RATIO Student/faculty ratio
  • DONATE Pct.alumni who donate
  • NEW10 Pct. new students from top 10% of H.S. class - % студентов из топ 10% своей старшей школы
  • NEW25 Pct. new students from top 25% of H.S. class - % студентов из топ 25% своей старшей школы

Качественные признаки:

  • FICE - Federal ID Number
  • …1 - Название университета
  • PPIND Public/private indicator (public=1, private=2)
  • STATE State (postal code)
  • TYPE - I (можно удалить)

Всего 49 признаков.

Если признаков очень много, то отобрать признаки (примерно 7-10)…

Отберём из условий задачи и по смыслу

  1. PPIND (для группировки)

Из условий задачи: 2. ADD_FEE - дополнительные выплаты 3. BOOK - плата за бронирование 4. NEW10 (зависимая переменная, также NEW25)

NEW10 зависит от AVRMATH, MATH_1, MATH_3, AVRVERB, VERB_1, VERB_3, AVR_ACT, ACT_1, ACT_3 и AVRCOMB. Признаков много, так что выберем самый обобщающий - AVRCOMB 5. AVRCOMB

  1. SF_RATIO

  2. PH_D

  3. GRADUAT 6-8 из предположения, что они взаимосвязаны

  4. R_B_COST

  5. IN_STATE

  6. OUT_STAT 9-11 - зависимость стоимости проживания (как окажется дальше из matrix plot, IN_STATE - фиктивный признак, то есть не информативный)

Определить вид признаков (колич., порядковые, качеств.)…

Проделаем это для рассматриваемых 11 признаков.

Setup

colNames<-c("ADD_FEE", "BOOK", "R_B_COST","IN_STATE","OUT_STAT","NEW10","AVRCOMB","SF_RATIO","PH_D","GRADUAT","PPIND")
colFin<-c("ADD_FEE", "BOOK", "PPIND")
colCost<-c("R_B_COST", "IN_STATE", "OUT_STAT", "PPIND")
colNew<-c("NEW10", "AVRCOMB", "PPIND")
colStud<-c("SF_RATIO", "PH_D", "GRADUAT", "PPIND")
col_I_comp<-col_I_sn[colNames]
col_I_comp_log<-col_I_sn_log[c(paste0("log-", colNames)[-length(colNames)], "PPIND")]
head(col_I_comp_log)|>rmarkdown::paged_table()
unique_info <- function(column) {
  column<-na.omit(column)
  list(ratio=length(unique(column)) / length(column), unique=length(unique(column)), total=length(column), moda = max(table(column)))
}

ratios <- sapply(col_I_comp, unique_info)
ratios
##        ADD_FEE BOOK      R_B_COST  IN_STATE  OUT_STAT  NEW10 AVRCOMB  
## ratio  0.875   0.3793103 0.9431818 0.9053254 0.9314286 0.425 0.7889908
## unique 119     66        166       153       163       68    86       
## total  136     174       176       169       175       160   109      
## moda   4       29        3         4         4         10    3        
##        SF_RATIO  PH_D      GRADUAT   PPIND     
## ratio  0.6363636 0.1976048 0.3918129 0.01136364
## unique 112       33        67        2         
## total  176       167       171       176       
## moda   5         10        7         123

Количественные признаки

Непрерывные признаки

  1. ADD_FEE
  2. R_B_COST
  3. IN_STATE
  4. OUT_STATE
  5. AVRCOMB
  6. SF_RATIO
  7. GRADUAT

Дискретные признаки

  1. BOOK (непрерывный, но дискретный из-за моды. Может быть округление)
  2. PH_D (считаем из-за плохой точности и округления данных)
  3. NEW10 (по ratio и moda, низкая точность)

Качественные признаки

  • PPIND - private/public.

Если признак порядковый и для него использованы текстовые метки, то проверить, что кодировка текстовых меток соответствует их естественному порядку.


Построить matrix plot (pairs plot), его долго разглядывать с точки зрения outliers, неоднородностей, вида распределений, вида зависимостей (линейные/нелинейные) и пр.

Первое построение pairs plot. Общие наблюдения.

ggpairs(col_I_comp, aes(alpha = 0.5, color=PPIND),
        diag = list(continuous = wrap("barDiag", bins = 20)),
        upper = list(continuous = wrap("cor", size = 2)))

Наблюдения:

  1. IN_STATE и OUT_STAT для private выстраивает линейную зависимость, для public - близкая к линейной зависимости. Отличительная особенность признака IN_STATE от OUT_STAT - наличие других факторов на образование цены обучения. OUT_STATE более приложим для анализа данных, поэтому IN_STATE можно считать фиктивным признаком.
  2. ADD_FEE похоже на log-normal. Заметны Outliers.
Outliers, ADD_FEE:
col_I_sn|>filter(ADD_FEE>3000)
## # A tibble: 7 × 49
##   ...1     PPIND  FICE STATE TYPE  AVRMATH AVRVERB AVRCOMB AVR_ACT MATH_1 MATH_3
##   <chr>    <chr> <dbl> <chr> <chr>   <dbl>   <dbl>   <dbl>   <dbl>  <dbl>  <dbl>
## 1 Univers… publ…  1313 CA    I         585     485    1070      24    520    650
## 2 Univers… publ…  1315 CA    I         621     520    1141      NA    558    700
## 3 Univers… publ…  1316 CA    I         542     440     982      23    470    620
## 4 Univers… publ…  1317 CA    I         616     515    1131      26    560    680
## 5 Univers… publ…  1320 CA    I         547     459    1006      NA    520    610
## 6 Univers… publ…  1321 CA    I         554     505    1059      NA    490    620
## 7 Univers… publ…  2221 MA    I          NA      NA      NA      NA    460    590
## # ℹ 38 more variables: VERB_1 <dbl>, VERB_3 <dbl>, ACT_1 <dbl>, ACT_3 <dbl>,
## #   APP_REC <dbl>, APP_ACC <dbl>, NEW_STUD <dbl>, NEW10 <dbl>, NEW25 <dbl>,
## #   FULLTIME <dbl>, PARTTIME <dbl>, IN_STATE <dbl>, OUT_STAT <dbl>,
## #   R_B_COST <dbl>, ROOM <dbl>, BOARD <dbl>, ADD_FEE <dbl>, BOOK <dbl>,
## #   PERSONAL <dbl>, PH_D <dbl>, TERM_D <dbl>, SF_RATIO <dbl>, DONATE <dbl>,
## #   INSTRUCT <dbl>, GRADUAT <dbl>, SAL_FULL <dbl>, SAL_AC <dbl>, SAL_AS <dbl>,
## #   SAL_ALL <dbl>, COMP_FUL <dbl>, COMP_AC <dbl>, COMP_AS <dbl>, …

University of California, Massachusetts - особые индивиды, в связи с высоким ВВП штатов.

  1. По histogram plot неоднородности (по PPIND) в R_B_COST, IN_STATE, OUT_STAT, SF_RATIO, GRADUAT.

Видим, что NEW10 зависит от PH_D, GRADUAT, AVRCOMB, отдельно для private зависит от OUT_STAT, IN_STATE, для public - от ADD_FEE, BOOK, R_B_COST

Если есть сильно несимметричные (с хвостом вправо) распределения на положительной полуоси, то…

Второе построение pairs plot.

Прологарифмируем ADD_FEE

col_I_mixed<-col_I_comp|>mutate(ADD_FEE=col_I_comp_log$`log-ADD_FEE`)|>rename(`log-ADD_FEE`=ADD_FEE)|>dplyr::select(-IN_STATE)
head(col_I_mixed)
## # A tibble: 6 × 10
##   `log-ADD_FEE`  BOOK R_B_COST OUT_STAT NEW10 AVRCOMB SF_RATIO  PH_D GRADUAT
##           <dbl> <dbl>    <dbl>    <dbl> <dbl>   <dbl>    <dbl> <dbl>   <dbl>
## 1         NA      600     3933     6300    25    1076     16.7    85      69
## 2         NA      580     3530     5424    NA      NA     17.3    80      50
## 3          5.67   750     5175     4440    24      NA      6.7    96      33
## 4          5.04   650     3590     5226    NA     961     10      67      NA
## 5          4.19   700     4850     7434    24     974     18.9    88      48
## 6          4.19   620     3728     6746    23     939     21.7    78      41
## # ℹ 1 more variable: PPIND <chr>
ggpairs(col_I_mixed, aes(alpha = 0.5, color=PPIND),
        diag = list(continuous = wrap("barDiag", bins = 20)),
        upper = list(continuous = wrap("cor", size = 2)))

Наблюдения:

  • Неоднородности (по PPIND) R_B_COST, OUT_STAT, NEW10.
  • Близкие к линейномй зависимости NEW10/ AVRCOMB.
  • Outliers, BOOK
Outliers, BOOK

Особые индивиды

col_I_sn|>filter(BOOK>=900 | BOOK <= 350)
## # A tibble: 4 × 49
##   ...1     PPIND  FICE STATE TYPE  AVRMATH AVRVERB AVRCOMB AVR_ACT MATH_1 MATH_3
##   <chr>    <chr> <dbl> <chr> <chr>   <dbl>   <dbl>   <dbl>   <dbl>  <dbl>  <dbl>
## 1 Howard … priv…  1448 DC    I         457     426     883      20    390    510
## 2 Tulane … priv…  2029 LA    I          NA      NA      NA      NA    560    670
## 3 Hofstra… priv…  2732 NY    I         540     470    1010      24    480    590
## 4 Renssel… priv…  2803 NY    I          NA      NA      NA      28    610    710
## # ℹ 38 more variables: VERB_1 <dbl>, VERB_3 <dbl>, ACT_1 <dbl>, ACT_3 <dbl>,
## #   APP_REC <dbl>, APP_ACC <dbl>, NEW_STUD <dbl>, NEW10 <dbl>, NEW25 <dbl>,
## #   FULLTIME <dbl>, PARTTIME <dbl>, IN_STATE <dbl>, OUT_STAT <dbl>,
## #   R_B_COST <dbl>, ROOM <dbl>, BOARD <dbl>, ADD_FEE <dbl>, BOOK <dbl>,
## #   PERSONAL <dbl>, PH_D <dbl>, TERM_D <dbl>, SF_RATIO <dbl>, DONATE <dbl>,
## #   INSTRUCT <dbl>, GRADUAT <dbl>, SAL_FULL <dbl>, SAL_AC <dbl>, SAL_AS <dbl>,
## #   SAL_ALL <dbl>, COMP_FUL <dbl>, COMP_AC <dbl>, COMP_AS <dbl>, …

Если есть outliers, то попробовать объяснить причину (ошибка в данных, особые индивиды) и удалить их.

BOOK

ggpairs(col_I_mixed[,c("log-ADD_FEE", "BOOK", "PPIND")], aes(alpha = 0.5, color=PPIND),
        diag = list(continuous = wrap("barDiag", bins = 20)),
        upper = list(continuous = wrap("cor", size = 2)))
## Warning: Removed 40 rows containing non-finite outside the scale range
## (`stat_bin()`).
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 40 rows containing missing values
## Warning: Removed 40 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
## Warning: Removed 40 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_bin()`).
## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 40 rows containing non-finite outside the scale range
## (`stat_bin()`).
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_bin()`).

col_I_sn|>filter(BOOK>=900 | BOOK <= 350)
## # A tibble: 4 × 49
##   ...1     PPIND  FICE STATE TYPE  AVRMATH AVRVERB AVRCOMB AVR_ACT MATH_1 MATH_3
##   <chr>    <chr> <dbl> <chr> <chr>   <dbl>   <dbl>   <dbl>   <dbl>  <dbl>  <dbl>
## 1 Howard … priv…  1448 DC    I         457     426     883      20    390    510
## 2 Tulane … priv…  2029 LA    I          NA      NA      NA      NA    560    670
## 3 Hofstra… priv…  2732 NY    I         540     470    1010      24    480    590
## 4 Renssel… priv…  2803 NY    I          NA      NA      NA      28    610    710
## # ℹ 38 more variables: VERB_1 <dbl>, VERB_3 <dbl>, ACT_1 <dbl>, ACT_3 <dbl>,
## #   APP_REC <dbl>, APP_ACC <dbl>, NEW_STUD <dbl>, NEW10 <dbl>, NEW25 <dbl>,
## #   FULLTIME <dbl>, PARTTIME <dbl>, IN_STATE <dbl>, OUT_STAT <dbl>,
## #   R_B_COST <dbl>, ROOM <dbl>, BOARD <dbl>, ADD_FEE <dbl>, BOOK <dbl>,
## #   PERSONAL <dbl>, PH_D <dbl>, TERM_D <dbl>, SF_RATIO <dbl>, DONATE <dbl>,
## #   INSTRUCT <dbl>, GRADUAT <dbl>, SAL_FULL <dbl>, SAL_AC <dbl>, SAL_AS <dbl>,
## #   SAL_ALL <dbl>, COMP_FUL <dbl>, COMP_AC <dbl>, COMP_AS <dbl>, …

Может быть как низкой, так и высокой в зависимости от штата, политики университета и скрытых факторов.

ADD_FEE

col_I_sn|>filter(ADD_FEE>3000)
## # A tibble: 7 × 49
##   ...1     PPIND  FICE STATE TYPE  AVRMATH AVRVERB AVRCOMB AVR_ACT MATH_1 MATH_3
##   <chr>    <chr> <dbl> <chr> <chr>   <dbl>   <dbl>   <dbl>   <dbl>  <dbl>  <dbl>
## 1 Univers… publ…  1313 CA    I         585     485    1070      24    520    650
## 2 Univers… publ…  1315 CA    I         621     520    1141      NA    558    700
## 3 Univers… publ…  1316 CA    I         542     440     982      23    470    620
## 4 Univers… publ…  1317 CA    I         616     515    1131      26    560    680
## 5 Univers… publ…  1320 CA    I         547     459    1006      NA    520    610
## 6 Univers… publ…  1321 CA    I         554     505    1059      NA    490    620
## 7 Univers… publ…  2221 MA    I          NA      NA      NA      NA    460    590
## # ℹ 38 more variables: VERB_1 <dbl>, VERB_3 <dbl>, ACT_1 <dbl>, ACT_3 <dbl>,
## #   APP_REC <dbl>, APP_ACC <dbl>, NEW_STUD <dbl>, NEW10 <dbl>, NEW25 <dbl>,
## #   FULLTIME <dbl>, PARTTIME <dbl>, IN_STATE <dbl>, OUT_STAT <dbl>,
## #   R_B_COST <dbl>, ROOM <dbl>, BOARD <dbl>, ADD_FEE <dbl>, BOOK <dbl>,
## #   PERSONAL <dbl>, PH_D <dbl>, TERM_D <dbl>, SF_RATIO <dbl>, DONATE <dbl>,
## #   INSTRUCT <dbl>, GRADUAT <dbl>, SAL_FULL <dbl>, SAL_AC <dbl>, SAL_AS <dbl>,
## #   SAL_ALL <dbl>, COMP_FUL <dbl>, COMP_AC <dbl>, COMP_AS <dbl>, …

University of California, Massachusetts - особые индивиды (из-за ВВП).

AVRCOMB/NEW10

ggpairs(col_I_mixed[colNew], aes(alpha = 0.5, color=PPIND),
        diag = list(continuous = wrap("barDiag", bins = 20)),
        upper = list(continuous = wrap("cor", size = 2)))
## Warning: Removed 16 rows containing non-finite outside the scale range
## (`stat_bin()`).
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 75 rows containing missing values
## Warning: Removed 16 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
## Warning: Removed 75 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 67 rows containing non-finite outside the scale range
## (`stat_bin()`).
## Warning: Removed 67 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 16 rows containing non-finite outside the scale range
## (`stat_bin()`).
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 67 rows containing non-finite outside the scale range
## (`stat_bin()`).

col_I_sn|>filter(AVRCOMB<=1200 & NEW10>=75)
## # A tibble: 8 × 49
##   ...1     PPIND  FICE STATE TYPE  AVRMATH AVRVERB AVRCOMB AVR_ACT MATH_1 MATH_3
##   <chr>    <chr> <dbl> <chr> <chr>   <dbl>   <dbl>   <dbl>   <dbl>  <dbl>  <dbl>
## 1 Univers… publ…  1313 CA    I         585     485    1070      24    520    650
## 2 Univers… publ…  1314 CA    I         574     455    1029      NA    510    640
## 3 Univers… publ…  1315 CA    I         621     520    1141      NA    558    700
## 4 Univers… publ…  1316 CA    I         542     440     982      23    470    620
## 5 Univers… publ…  1317 CA    I         616     515    1131      26    560    680
## 6 Univers… publ…  1320 CA    I         547     459    1006      NA    520    610
## 7 Univers… publ…  1321 CA    I         554     505    1059      NA    490    620
## 8 Univers… publ…  2974 NC    I         594     527    1121      NA    530    530
## # ℹ 38 more variables: VERB_1 <dbl>, VERB_3 <dbl>, ACT_1 <dbl>, ACT_3 <dbl>,
## #   APP_REC <dbl>, APP_ACC <dbl>, NEW_STUD <dbl>, NEW10 <dbl>, NEW25 <dbl>,
## #   FULLTIME <dbl>, PARTTIME <dbl>, IN_STATE <dbl>, OUT_STAT <dbl>,
## #   R_B_COST <dbl>, ROOM <dbl>, BOARD <dbl>, ADD_FEE <dbl>, BOOK <dbl>,
## #   PERSONAL <dbl>, PH_D <dbl>, TERM_D <dbl>, SF_RATIO <dbl>, DONATE <dbl>,
## #   INSTRUCT <dbl>, GRADUAT <dbl>, SAL_FULL <dbl>, SAL_AC <dbl>, SAL_AS <dbl>,
## #   SAL_ALL <dbl>, COMP_FUL <dbl>, COMP_AC <dbl>, COMP_AS <dbl>, …

University of California, North Carolina. Низкий порог для поступления?

PH_D

ggpairs(col_I_mixed[colStud], aes(alpha = 0.5, color=PPIND),
        diag = list(continuous = wrap("barDiag", bins = 20)),
        upper = list(continuous = wrap("cor", size = 2)))
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 9 rows containing missing values
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 5 rows containing missing values
## Warning: Removed 9 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 9 rows containing non-finite outside the scale range
## (`stat_bin()`).
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removed 14 rows containing missing values
## Warning: Removed 9 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
## Warning: Removed 5 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 14 rows containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 5 rows containing non-finite outside the scale range
## (`stat_bin()`).
## Warning: Removed 5 rows containing non-finite outside the scale range
## (`stat_boxplot()`).
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 9 rows containing non-finite outside the scale range
## (`stat_bin()`).
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 5 rows containing non-finite outside the scale range
## (`stat_bin()`).

col_I_sn|>filter(PH_D<=70)
## # A tibble: 5 × 49
##   ...1     PPIND  FICE STATE TYPE  AVRMATH AVRVERB AVRCOMB AVR_ACT MATH_1 MATH_3
##   <chr>    <chr> <dbl> <chr> <chr>   <dbl>   <dbl>   <dbl>   <dbl>  <dbl>  <dbl>
## 1 Univers… publ…  1063 AK    I         499     462     961      22     NA     NA
## 2 Ball St… publ…  1786 IN    I         473     422     895      21    400    530
## 3 Univers… publ…  2031 LA    I          NA      NA      NA      19     NA     NA
## 4 Texas S… publ…  3642 TX    I         390     420     810      19     NA     NA
## 5 Texas W… publ…  3646 TX    I          NA      NA      NA      NA    350    480
## # ℹ 38 more variables: VERB_1 <dbl>, VERB_3 <dbl>, ACT_1 <dbl>, ACT_3 <dbl>,
## #   APP_REC <dbl>, APP_ACC <dbl>, NEW_STUD <dbl>, NEW10 <dbl>, NEW25 <dbl>,
## #   FULLTIME <dbl>, PARTTIME <dbl>, IN_STATE <dbl>, OUT_STAT <dbl>,
## #   R_B_COST <dbl>, ROOM <dbl>, BOARD <dbl>, ADD_FEE <dbl>, BOOK <dbl>,
## #   PERSONAL <dbl>, PH_D <dbl>, TERM_D <dbl>, SF_RATIO <dbl>, DONATE <dbl>,
## #   INSTRUCT <dbl>, GRADUAT <dbl>, SAL_FULL <dbl>, SAL_AC <dbl>, SAL_AS <dbl>,
## #   SAL_ALL <dbl>, COMP_FUL <dbl>, COMP_AC <dbl>, COMP_AS <dbl>, …

Незначительные выбросы - особые индивиды. Могут быть связаны с низким GRADUAT и высоким SF_RATIO.

OUT_STAT

ggpairs(col_I_mixed[,c("R_B_COST", "OUT_STAT", "PPIND")], aes(alpha = 0.5, color=PPIND),
        diag = list(continuous = wrap("barDiag", bins = 20)),
        upper = list(continuous = wrap("cor", size = 2)))
## Warning in ggally_statistic(data = data, mapping = mapping, na.rm = na.rm, :
## Removing 1 row that contained a missing value
## Warning: Removed 1 row containing missing values or values outside the scale range
## (`geom_point()`).
## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).
## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_boxplot()`).
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_bin()`).

col_I_sn|>filter(OUT_STAT>=14000 & PPIND=="public")
## # A tibble: 2 × 49
##   ...1     PPIND  FICE STATE TYPE  AVRMATH AVRVERB AVRCOMB AVR_ACT MATH_1 MATH_3
##   <chr>    <chr> <dbl> <chr> <chr>   <dbl>   <dbl>   <dbl>   <dbl>  <dbl>  <dbl>
## 1 Univers… publ…  3696 VT    I         553     488    1041      NA    500    610
## 2 Univers… publ…  9092 MI    I         634     543    1177      27    580    700
## # ℹ 38 more variables: VERB_1 <dbl>, VERB_3 <dbl>, ACT_1 <dbl>, ACT_3 <dbl>,
## #   APP_REC <dbl>, APP_ACC <dbl>, NEW_STUD <dbl>, NEW10 <dbl>, NEW25 <dbl>,
## #   FULLTIME <dbl>, PARTTIME <dbl>, IN_STATE <dbl>, OUT_STAT <dbl>,
## #   R_B_COST <dbl>, ROOM <dbl>, BOARD <dbl>, ADD_FEE <dbl>, BOOK <dbl>,
## #   PERSONAL <dbl>, PH_D <dbl>, TERM_D <dbl>, SF_RATIO <dbl>, DONATE <dbl>,
## #   INSTRUCT <dbl>, GRADUAT <dbl>, SAL_FULL <dbl>, SAL_AC <dbl>, SAL_AS <dbl>,
## #   SAL_ALL <dbl>, COMP_FUL <dbl>, COMP_AC <dbl>, COMP_AS <dbl>, …
col_I_sn|>filter(OUT_STAT<=10500 & PPIND=="private")
## # A tibble: 6 × 49
##   ...1     PPIND  FICE STATE TYPE  AVRMATH AVRVERB AVRCOMB AVR_ACT MATH_1 MATH_3
##   <chr>    <chr> <dbl> <chr> <chr>   <dbl>   <dbl>   <dbl>   <dbl>  <dbl>  <dbl>
## 1 Univers… priv…  1431 DE    I          NA      NA    1035      NA    500    610
## 2 Howard … priv…  1448 DC    I         457     426     883      20    390    510
## 3 Rice Un… priv…  3604 TX    I          NA      NA      NA      NA    665    763
## 4 Texas C… priv…  3636 TX    I          NA      NA      NA      NA    460    600
## 5 Brigham… priv…  3670 UT    I          NA      NA      NA      NA     NA     NA
## 6 Baylor … priv…  6967 TX    I          NA      NA      NA      NA    500    620
## # ℹ 38 more variables: VERB_1 <dbl>, VERB_3 <dbl>, ACT_1 <dbl>, ACT_3 <dbl>,
## #   APP_REC <dbl>, APP_ACC <dbl>, NEW_STUD <dbl>, NEW10 <dbl>, NEW25 <dbl>,
## #   FULLTIME <dbl>, PARTTIME <dbl>, IN_STATE <dbl>, OUT_STAT <dbl>,
## #   R_B_COST <dbl>, ROOM <dbl>, BOARD <dbl>, ADD_FEE <dbl>, BOOK <dbl>,
## #   PERSONAL <dbl>, PH_D <dbl>, TERM_D <dbl>, SF_RATIO <dbl>, DONATE <dbl>,
## #   INSTRUCT <dbl>, GRADUAT <dbl>, SAL_FULL <dbl>, SAL_AC <dbl>, SAL_AS <dbl>,
## #   SAL_ALL <dbl>, COMP_FUL <dbl>, COMP_AC <dbl>, COMP_AS <dbl>, …

Цены на обучение варьируются по штатам и самим университетам.


Если есть неоднородности (например, видны два облака точек), то объяснить причину (найти категоризующую переменную, объясняющую эту неоднородность).

ggpairs(col_I_mixed, aes(alpha = 0.5, color=PPIND),
        diag = list(continuous = wrap("barDiag", bins = 20)),
        upper = list(continuous = wrap("cor", size = 2)))

Выделим неоднородности R_B_COST, OUT_STAT. Явное разделение на два облака точек для этих признаков может быть связано с разницей финансирования. Для public - государство, штат и граждане, для private - только граждане.


В дальнейшем вид pairs plots, распределения признаков и корреляции анализировать отдельно для неоднородных групп.


Посмотрите также на descriptive statistics с точки зрения минимумов-максимумов, асимметрии, эксцесса и пр.

results <- describeBy(col_I_mixed, group=col_I_comp$PPIND)
df_results <- map_dfr(results, ~as.data.frame(.x), .id = "group")
kable(df_results, "html") %>%
  kable_styling(full_width = FALSE) %>%
  scroll_box(width = "100%", height = "400px")
group vars n mean sd median trimmed mad min max range skew kurtosis se
log-ADD_FEE…1 private 1 46 5.838785 0.7936613 6.003811 5.872662 0.7083206 3.688880 7.515345 3.826465 -0.4493357 0.0057261 0.1170190
BOOK…2 private 2 53 610.075472 158.1109456 600.000000 597.534884 148.2600000 300.000000 1230.000000 930.000000 1.2325851 3.0965764 21.7182087
R_B_COST…3 private 3 53 5839.811321 997.0704793 5975.000000 5922.255814 816.9126000 3320.000000 7425.000000 4105.000000 -0.7233444 -0.2295138 136.9581633
OUT_STAT…4 private 4 53 15747.358491 3924.9053349 17020.000000 16306.186047 2772.4620000 2340.000000 20100.000000 17760.000000 -1.3082365 1.3192370 539.1272103
NEW10…5 private 5 52 57.788461 24.3796255 58.000000 57.761905 33.3585000 16.000000 98.000000 82.000000 0.0381849 -1.4295240 3.3808458
AVRCOMB…6 private 6 29 1177.034483 144.2860459 1199.000000 1177.680000 161.6034000 883.000000 1410.000000 527.000000 -0.0893581 -1.2296542 26.7932461
SF_RATIO…7 private 7 53 9.813208 4.3268854 9.200000 9.555814 4.8925800 2.900000 20.500000 17.600000 0.4758614 -0.7282486 0.5943434
PH_D…8 private 8 49 89.591837 7.8658502 91.000000 90.317073 7.4130000 71.000000 99.000000 28.000000 -0.8078233 -0.4898133 1.1236929
GRADUAT…9 private 9 52 78.711539 15.7534142 78.500000 80.023809 17.7912000 33.000000 99.000000 66.000000 -0.6127581 -0.1790876 2.1846055
PPIND…10 private 10 53 1.000000 0.0000000 1.000000 1.000000 0.0000000 1.000000 1.000000 0.000000 NaN NaN 0.0000000
log-ADD_FEE…11 public 1 90 6.014273 1.1296149 6.093567 6.024473 0.8723448 2.995732 8.383433 5.387701 -0.1980167 0.2896088 0.1190719
BOOK…12 public 2 121 600.008264 101.1621220 600.000000 591.103093 111.1950000 400.000000 858.000000 458.000000 0.6447669 -0.2889320 9.1965565
R_B_COST…13 public 3 123 3977.073171 842.4904744 3811.000000 3911.979798 722.0262000 2082.000000 6607.000000 4525.000000 0.7886999 0.7902446 75.9648078
OUT_STAT…14 public 4 122 7657.549180 2252.3124723 7446.500000 7517.163265 1939.2408000 2279.000000 15732.000000 13453.000000 0.7771996 1.2996173 203.9147900
NEW10…15 public 5 108 33.620370 21.8669373 26.000000 29.920455 11.8608000 8.000000 95.000000 87.000000 1.5420644 1.5799521 2.1041470
AVRCOMB…16 public 6 80 1014.662500 84.6586510 997.000000 1010.687500 75.6126000 810.000000 1240.000000 430.000000 0.4180047 0.1997793 9.4651249
SF_RATIO…17 public 7 123 16.129268 3.8302021 16.500000 16.178788 4.1512800 6.700000 24.700000 18.000000 -0.0994206 -0.5699091 0.3453577
PH_D…18 public 8 118 84.110169 7.7985139 85.000000 84.406250 7.4130000 63.000000 99.000000 36.000000 -0.3698537 -0.3490528 0.7179114
GRADUAT…19 public 9 119 54.731092 15.5730632 54.000000 54.412371 16.3086000 10.000000 95.000000 85.000000 0.1127780 -0.1892703 1.4275804
PPIND…20 public 10 123 2.000000 0.0000000 2.000000 2.000000 0.0000000 2.000000 2.000000 0.000000 NaN NaN 0.0000000

Наблюдения:

Mean, median, skew:

NEW10 и AVRCOMB имеют skew, близкий к 0 (private), SF_RATIO, GRADUAT, log-ADD_FEE околонулевой skew (public).

Kurtosis:

Близкий к 0 (private): log-ADD_FEE (но skew ~ -0.44).

2. О виде распределений и о сравнении распределений

Первые два пункта индивидуального задания нужно делать не по указанному порядку, а как того требует логика статистики. Чтобы сравнивать выборки по t-критерию, нужно знать о том, близки ли распределения в сравниваемых группах к нормальным или хотя бы к унимодальным и симметричным. Чтобы проверять распределения признаков на нормальность, нужно знать, что рассматривается однородная выборка.

col_split_mixed<-split(col_I_mixed, col_I_mixed$PPIND)
names(col_split_mixed)<-c("private", "public")
colMixed<-names(col_I_mixed)
colMixed<-colMixed[-length(colMixed)]
colMixed
## [1] "log-ADD_FEE" "BOOK"        "R_B_COST"    "OUT_STAT"    "NEW10"      
## [6] "AVRCOMB"     "SF_RATIO"    "PH_D"        "GRADUAT"

Проверка на нормальность по хи-квадрат критерию

chi_squared_normality_test <- function(data, min_expected = 4.9) {
  data <- na.omit(data)
  n <- length(data)
  if (n < 10) stop("Sample size too small for this test. Need at least 10 observations.")

  mu <- mean(data)
  sigma <- sd(data)

   # Calculate quantiles for breaks
   probs <- seq(5/n, 1 - 5/n, by = 5/n)
   breaks <- qnorm(probs, mean = mu, sd = sigma)

   # Add extreme quantiles
   breaks <- c(qnorm(0), breaks, qnorm(1))

  #Observed Frequencies
  observed <- hist(data, breaks = breaks, plot = FALSE)$counts

  #Expected Frequencies (using pnorm)
    expected <- diff(pnorm(breaks, mean = mu, sd = sigma)) * n

  chi_squared <- sum((observed - expected)^2 / expected)
  df <- length(expected) - 3  #Corrected degrees of freedom

  p_value <- pchisq(chi_squared, df = df, lower.tail = FALSE)

  return(list(p_value = p_value, statistic = chi_squared, df = df))
}

Проверка на достоверность (p-value должно иметь равномерное распределение на (0, 1))

n = 500
test<-replicate(5000, chi_squared_normality_test(rnorm(n))$p_value)
quantile(test)
##           0%          25%          50%          75%         100% 
## 0.0002821321 0.2499259663 0.5038633005 0.7519041956 0.9996275164

Тесты на нормальность

perform_normality_tests <- function(df, cols) {
  # Input validation (remains the same)
    # Function to perform a single normality test (simplified)
  run_test <- function(data, test_name) {
    if (test_name == "Chi-squared") {
      chisq_test<-chi_squared_normality_test(data)
      return(chisq_test[c("p_value", "statistic", "df")])
    }
    test_result <- switch(test_name,
                          "Lilliefors" = lillie.test(data),
                          "Anderson-Darling" = ad.test(data),
                          "Shapiro-Wilk" = shapiro.test(data)
    )
    return(list(p_value = test_result$p.value, statistic = test_result$statistic, df=length(na.omit(data))))  # Removed test_name
  }
     results <- lapply(cols, function(col) {
    data <- df[[col]]
    n_na <- sum(is.na(data))
    if (n_na > 0) {
      warning(paste0("Removed ", n_na, " NA values from column '", col, "' before testing."))
      data <- na.omit(data)
    }
    test_results <- lapply(c("Lilliefors", "Anderson-Darling", "Shapiro-Wilk", "Chi-squared"), function(test) run_test(data, test))
    list(column = col, tests = test_results)
  })

  results_df <- do.call(rbind, lapply(results, function(x) {
    test_names <- c("Lilliefors", "Anderson-Darling", "Shapiro-Wilk", "Chi-squared")
    p_values <- sapply(x$tests, "[[", "p_value")
    data.frame(
      Column = x$column,
      Test = test_names, #Directly use test_names vector
      Statistic = sapply(x$tests, "[[", "statistic"),
      P_value = p_values,
      Size = sapply(x$tests, "[[", "df"),
      Significance = sapply(p_values, function(p) ifelse(p < 0.05, "\u2714", "\u2716")),
      stringsAsFactors = FALSE
    )
  }))

  plots <- lapply(cols, function(col) {
    data <- df[[col]]
    data <- na.omit(data)
    ggqqplot(data, title = paste("Normal Probability Plot for", col))

  })
  names(plots)<-cols
  list(results = results_df, plots = plots)
}

Так как визуально однородность при предварительном анализе была уже исследована, то можно начинать с анализа вида распределения признаков, возможно, по группам. Сюда входит: normal probability plot (что это такое?), проверка по критериям Лиллиефорса, AD, хи-квадрат, Шапиро-Уилка. По критерию хи-квадрат, а также визуально по PP-plot можно проверить и гипотезы о согласии с другими распределениями, например, логнормальным. Задаваемые вопросы: чем отличается критерий Лиллиефорса от критерия Колмогорова, в чем отличие AD, как выглядят статистики критериев, что такое PP-plot и normal probability plot, почему естественно при рисовании normal probability plot одновременно смотреть на результаты критерия Шапиро-Уилка.

Public: проверка на нормальность, тесты

# colSums(is.na(col_split_mixed$public))
col_norm_public<-col_split_mixed$public|>perform_normality_tests(colMixed)
## Warning in FUN(X[[i]], ...): Removed 33 NA values from column 'log-ADD_FEE'
## before testing.
## Warning in FUN(X[[i]], ...): Removed 2 NA values from column 'BOOK' before
## testing.
## Warning in FUN(X[[i]], ...): Removed 1 NA values from column 'OUT_STAT' before
## testing.
## Warning in FUN(X[[i]], ...): Removed 15 NA values from column 'NEW10' before
## testing.
## Warning in FUN(X[[i]], ...): Removed 43 NA values from column 'AVRCOMB' before
## testing.
## Warning in FUN(X[[i]], ...): Removed 5 NA values from column 'PH_D' before
## testing.
## Warning in FUN(X[[i]], ...): Removed 4 NA values from column 'GRADUAT' before
## testing.
col_norm_public$plots
## $`log-ADD_FEE`

## 
## $BOOK

## 
## $R_B_COST

## 
## $OUT_STAT

## 
## $NEW10

## 
## $AVRCOMB

## 
## $SF_RATIO

## 
## $PH_D

## 
## $GRADUAT

col_norm_public$results|>dplyr::select(-Statistic)|>rmarkdown::paged_table()

Наблюдения:

Хи-квадрат: не отвергаем AVRCOMB, SF_RATIO, PH_D, GRADUAT, R_B_COST, OUT_STAT,

Можно заметить “хвосты” на некоторых графиках.

Private: проверка на нормальность, тесты

col_norm_private<-col_split_mixed$private|>perform_normality_tests(colMixed)
## Warning in FUN(X[[i]], ...): Removed 7 NA values from column 'log-ADD_FEE'
## before testing.
## Warning in FUN(X[[i]], ...): Removed 1 NA values from column 'NEW10' before
## testing.
## Warning in FUN(X[[i]], ...): Removed 24 NA values from column 'AVRCOMB' before
## testing.
## Warning in FUN(X[[i]], ...): Removed 4 NA values from column 'PH_D' before
## testing.
## Warning in FUN(X[[i]], ...): Removed 1 NA values from column 'GRADUAT' before
## testing.
col_norm_private$plots
## $`log-ADD_FEE`

## 
## $BOOK

## 
## $R_B_COST

## 
## $OUT_STAT

## 
## $NEW10

## 
## $AVRCOMB

## 
## $SF_RATIO

## 
## $PH_D

## 
## $GRADUAT

col_norm_private$results|>dplyr::select(-Statistic)|>rmarkdown::paged_table()

Наблюдения:

Хи-квадрат: не отвергаем log-ADD_FEE, R_B_COST, NEW10, AVRCOMB, SF_RATIO, GRADUAT

Сначала имеет смысл посмотреть на сравнение распределений в группах с помощью ящиков с усами. С помощью ящиков с усами там, где групп больше двух, можно выбрать две из них, которые интересно сравнить с помощью критериев.

Совместно с 4.

Если в задании есть сравнение независимых выборок (точнее, распределений независимых случайных величин), то начинать нужно с t-критерия, который мощный против альтернатив, заключающихся в наиболее легко интерпретируемом сдвиге (т.е. разнице средних). Нужно не забыть, что у критерия есть варианты для модели с одинаковыми дисперсиями (получается точное p-value, которое может быть неправильным, если на самом деле дисперсии одинаковые) и с произвольными дисперсиями. В результате, получатся два критерия для гипотезы о равенстве средних и два критерия о равенстве разбросов. Нужно уметь объяснять, что это за критерии и при каких условиях их можно применять. Не забудьте, что при использовании асимптотических критериев нужно обращать внимание на объемы выборок. Сделайте выводы о том, для каких признаков есть разница в сдвиге.

Тесты на равенство дисперсий/матожиданий, тест Wilcoxon.

perform_tests <- function(data,
                          quantitative_columns,
                          categorical_column,
                          alpha = 0.05) {
  results <- list()
  
  for (col in quantitative_columns) {
    levene_test<-leveneTest(data[[col]] ~ data[[categorical_column]])
    var_test <- var.test(data[[col]] ~ data[[categorical_column]])
    wilcox_test <- wilcox.test(data[[col]]~data[[categorical_column]])
    t_test_equal <- t.test(data[[col]] ~ data[[categorical_column]], var.equal = TRUE)
    t_test_unequal <- t.test(data[[col]] ~ data[[categorical_column]], var.equal = FALSE)
    
    results[[col]] <- list(
      t_test_unequal_p_value = t_test_unequal$p.value,
      t_test_equal_p_value = t_test_equal$p.value,
      levene_test_p_value = levene_test$`Pr(>F)`[1],
      var_test_p_value = var_test$p.value,
      wilcox_test_p_value = wilcox_test$p.value
    )
  }
  
  return(results |> enframe(name = "attribute", value = "p_values") |>
  unnest_wider(p_values))
}
df_tests = perform_tests(col_I_mixed, colMixed, "PPIND")
## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.
## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.
## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.
## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.
## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.
## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.
## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.
## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.
## Warning in leveneTest.default(y = y, group = group, ...): group coerced to
## factor.
df_long <- col_I_mixed|>
  pivot_longer(cols = -PPIND, names_to = "attribute", values_to = "value")
# Объединяем результаты теста с основными данными
df_with_tests <- df_long %>%
  left_join(df_tests, by = "attribute")
df_with_tests_cost <- df_with_tests[df_with_tests$attribute%in%colCost[-length(colCost)],]
df_with_tests_new <- df_with_tests[df_with_tests$attribute%in%colNew[-length(colNew)],]
df_with_tests_stud <- df_with_tests[df_with_tests$attribute%in%colStud[-length(colStud)],]
df_with_tests_fin <- df_with_tests[df_with_tests$attribute%in%c("log-ADD_FEE", "BOOK"),]

p<-ggplot(df_with_tests_cost,
       aes(x = attribute,
           y = value,
           fill = PPIND)) +
  geom_boxplot() +
  labs(title = "Boxplot",
       x = "Attribute",
       y = "Value")+
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1), plot.caption = element_text(size = 8, hjust = 0)) +
  scale_fill_manual(values = c("private" = "red", "public" = "lightblue"))

AVRCOMB, NEW10.

p$data<-df_with_tests_new|>filter(attribute%in%"AVRCOMB")
p
## Warning: Removed 67 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

p$data<-df_with_tests_new|>filter(attribute%in%"NEW10")
p
## Warning: Removed 16 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

df_with_tests|>dplyr::filter(attribute%in%colNew[-length(colNew)])|>dplyr::select(-PPIND, -value)|>slice_head(n=2)|>rmarkdown::paged_table()

Наблюдения:

Можем сказать, что public и private имеют разные матожидания.

BOOK, log-ADD_FEE.

p$data<-df_with_tests_fin|>filter(attribute%in%"BOOK")
p
## Warning: Removed 2 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

p$data<-df_with_tests_fin|>filter(attribute%in%"log-ADD_FEE")
p
## Warning: Removed 40 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

df_with_tests|>dplyr::filter(attribute%in%c("log-ADD_FEE", "BOOK"))|>dplyr::select(-PPIND, -value)|>head(length(colFin)-1)|>rmarkdown::paged_table()

Наблюдения:

Исходя из p_value не отвергаем, что public и private имеют равные матожидания. Тест Вилкоксона также характеризует их однородность.

Задача 1.1 из .tsk

Проверим log-ADD_FEE и BOOK друг с другом

compare_column_tests <- function(column_1,
                          column_2,
                          alpha = 0.05) {
  levene_test<-leveneTest(column_2, column_1)
  var_test <- var.test(column_1, column_2)
  wilcox_test <- wilcox.test(column_1, column_2)
  t_test_equal <- t.test(column_1, column_2, var.equal = TRUE)
  t_test_unequal <- t.test(column_1, column_2, var.equal = FALSE)
  
  results <- data.frame(
    t_test_unequal_p_value = t_test_unequal$p.value,
    t_test_equal_p_value = t_test_equal$p.value,
    levene_test_p_value = levene_test$`Pr(>F)`[1],
    var_test_p_value = var_test$p.value,
    wilcox_test_p_value = wilcox_test$p.value
  )
}
af_pr<-col_I_sn$ADD_FEE[col_I_sn$PPIND=="private"]
af_pu<-col_I_sn$ADD_FEE[col_I_sn$PPIND=="public"]
b_pr<-col_split_mixed$private$BOOK
b_pu<-col_split_mixed$public$BOOK

print(compare_column_tests(af_pr, b_pr))
## Warning in leveneTest.default(column_2, column_1): column_1 coerced to factor.
## Warning in anova.lm(lm(resp ~ group)): применение F-критерия для в целом
## хорошей подгонки бессмысленно
##   t_test_unequal_p_value t_test_equal_p_value levene_test_p_value
## 1            0.005903872          0.003538701        2.925758e-75
##   var_test_p_value wilcox_test_p_value
## 1      1.42501e-07        1.678377e-05
print(compare_column_tests(af_pu, b_pu))
## Warning in leveneTest.default(column_2, column_1): column_1 coerced to factor.
##   t_test_unequal_p_value t_test_equal_p_value levene_test_p_value
## 1              0.1706573            0.1115642           0.9495775
##   var_test_p_value wilcox_test_p_value
## 1                0        0.0006091651
Наблюдения:

В общем случае тесты показывают, что разница есть в случае общественных вузов, но для частных вузов не отвергаем равенство дисперсий (из levene’s test) и матожиданий.

GRADUAT, PH_D, SF_RATIO

p$data<-df_with_tests_stud
p
## Warning: Removed 14 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

df_with_tests|>dplyr::filter(attribute%in%colStud[-length(colStud)])|>dplyr::select(-PPIND, -value)|>head(length(colStud)-1)|>rmarkdown::paged_table()

Наблюдения:

Не отвергаем равенство дисперсий. В случае с NEW10 по Levene Test отвергаем гипотезу, а по Var Test не отвергаем. Визуально, из-за выбросов Var Test становится мощнее, поэтому следует отвергнуть равенство дисперсий NEW10.

OUT_STAT, R_B_COST

p$data<-df_with_tests_cost
p
## Warning: Removed 1 row containing non-finite outside the scale range
## (`stat_boxplot()`).

df_with_tests|>dplyr::filter(attribute%in%colCost[-length(colCost)])|>dplyr::select(-PPIND, -value)|>head(length(colCost)-1)|>rmarkdown::paged_table()

Наблюдения:

Для R_B_COST не отвергается равенство дисперсий.

Далее, объясняете, в каких случаях (распределение далеко от нормального, могут быть выделяющиеся наблюдения, нас интересует другая характеристика положения) t-критерий не удовлетворителен и нужно переходить к непараметрическим критериям. Рассказываете, какой из непараметрических критериев является аналогом t-критерия, как он строится и против какой альтернативы мощный. Вы уже догадались, что это критерий Манна-Уитни, он же критерий Вилкоксона.

t-критерий неудовлетворителен в нескольких ситуациях:

  1. Малый объём выборки. Если форма распределения признака неизвестна или не нормальна, то t-критерий асимптотически точный.
  2. Наличие выбросов. t-критерий неусточив по отношению к выбросам.
  3. Другие цели эксперимента. t-критерий мощный против альтернатив о неравенстве матожиданий, его нельзя применить для построения других гипотез, например, гипотезы о равенстве дисперсий.

В качестве замены t-критерия можно воспользоваться непараметрическим критерием Манна-Уитни (Вилкоксона).

Смотрите на результаты применения критерия Манна-Уитни, сравниваете с результатами применения t-критерия. Проводите сравнительный анализ критериев с теоретической точки зрения (чем один лучше другого и чем хуже).

df_tests|>rmarkdown::paged_table()

Теоретическое сравнение:

Как непараметрический критерий, тест Вилкоксона уступает t-тесту в мощности против альтернативы о неравенстве матожиданий. С другой стороны, тест Вилкоксона устойчив к выбросам, мощный против альтернативы о неравенстве распределения двух случайных величин. Если известно, что две случайные величины имеют одинаковые формы распределения, то тест Вилкоксона автоматически ставит гипотезу о равенстве матожиданий.

Далее, переходите к критериям, которые умеют сравнивать не только характеристики положения, но и формы распределений (это критерий Колмогорова-Смирнова, например). Для каждого критерия (включая критерий Манна-Уитни), нужно уметь объяснять, что означают столбцы в таблицах результатов критериев. Также, при разных результатах проверки гипотезы о равенстве распределений нужно объяснять, почему один критерий, например, не отверг гипотезу, а другой – отверг.

В пунктах 1, 2 проверяли на нормальность распределения. Теперь проверим распределения по критерию Колмогорова-Смирнова.

Критерий Колмогорова-Смирнова (сравнение public и private)

numeric_cols <- col_I_mixed %>%
  dplyr::select(where(is.numeric)) %>%
  colnames()
results <- col_I_mixed %>%
    dplyr::select(PPIND, all_of(numeric_cols)) %>%
    pivot_longer(cols = !PPIND,
                     names_to = "attribute",
                    values_to = "value") %>%
    group_by(attribute) %>%
    summarize(
        test_result = list(
            ks.test(
                value[PPIND == "public"],
                value[PPIND == "private"]
            )
        ),
      .groups = 'drop'
    ) %>%
     mutate(
            p_value = map_dbl(test_result, ~ .x$p.value),
            statistic = map_dbl(test_result, ~ .x$statistic)
    ) %>%
    dplyr::select(-test_result)
results|>rmarkdown::paged_table()

Наблюдения:

  • Тест Колмогорова-Смирнова на равенство распределений public и private:

Не отвергаем:

  1. BOOK
  2. log-ADD_FEE

Об анализе зависимостей

Вспомните, какие бывают виды зависимостей и чем они измеряются, по каким формулам. Посмотрите на основе pairs plot, какие зависимости у вас в данных. Не забудьте, что при неоднородных данных изучать зависимости имеет смысл только внутри групп по-отдельности.

col_split_mixed<-split(col_I_mixed, col_I_mixed$PPIND)
names(col_split_mixed)<-c("private", "public")
col_split_mixed<-lapply(col_split_mixed, function(x){
  dplyr::select(x, -PPIND)
})
colMixed<-names(col_I_mixed)
colMixed<-colMixed[-length(colMixed)]
col_split_mixed$private|>rmarkdown::paged_table()
col_split_mixed$public|>rmarkdown::paged_table()

Private:

ggpairs(col_split_mixed$private, aes(alpha = 0.5),
        diag = list(continuous = wrap("barDiag", bins = 20)),
        upper = list(continuous = wrap("cor", size = 2)))

Наблюдения:

Можно заметить, что NEW10/AVRCOMB разбивается на две части (поскольку AVRCOMB имеет много пропусков и сам AVRCOMB подразумевает разделение на подгруппы).

col_I_sn[col_I_sn$PPIND=="private",]|>filter(NEW10<=50 & !is.na(AVRCOMB))|>rmarkdown::paged_table()
col_I_sn[col_I_sn$PPIND=="private",]|>filter(NEW10>50 & !is.na(AVRCOMB))|>rmarkdown::paged_table()

Среди частных университетов только один “Delaware University” имеет расхождения между IN_STATE и OUT_STAT. Сейчас идёт такое ценообразование.

University of Delaware.

University of Delaware’s tuition is $15,410 for in-state and $37,930 for out-of-state.

col_I_sn[col_I_sn$IN_STATE!=col_I_sn$OUT_STAT & col_I_sn$PPIND=="private",]
## # A tibble: 1 × 49
##   ...1     PPIND  FICE STATE TYPE  AVRMATH AVRVERB AVRCOMB AVR_ACT MATH_1 MATH_3
##   <chr>    <chr> <dbl> <chr> <chr>   <dbl>   <dbl>   <dbl>   <dbl>  <dbl>  <dbl>
## 1 Univers… priv…  1431 DE    I          NA      NA    1035      NA    500    610
## # ℹ 38 more variables: VERB_1 <dbl>, VERB_3 <dbl>, ACT_1 <dbl>, ACT_3 <dbl>,
## #   APP_REC <dbl>, APP_ACC <dbl>, NEW_STUD <dbl>, NEW10 <dbl>, NEW25 <dbl>,
## #   FULLTIME <dbl>, PARTTIME <dbl>, IN_STATE <dbl>, OUT_STAT <dbl>,
## #   R_B_COST <dbl>, ROOM <dbl>, BOARD <dbl>, ADD_FEE <dbl>, BOOK <dbl>,
## #   PERSONAL <dbl>, PH_D <dbl>, TERM_D <dbl>, SF_RATIO <dbl>, DONATE <dbl>,
## #   INSTRUCT <dbl>, GRADUAT <dbl>, SAL_FULL <dbl>, SAL_AC <dbl>, SAL_AS <dbl>,
## #   SAL_ALL <dbl>, COMP_FUL <dbl>, COMP_AC <dbl>, COMP_AS <dbl>, …

Public:

ggpairs(col_split_mixed$public, aes(alpha = 0.5),
        diag = list(continuous = wrap("barDiag", bins = 20)),
        upper = list(continuous = wrap("cor", size = 2)))

Наблюдения

На pairs plot можно заметить выбросы в признаках log-ADD_FEE, OUT_STAT, NEW10, которые влияют на коэффициент корреляции Пирсона.

Начинать нужно с анализа линейных зависимостей. На основе коэффициента корреляции Пирсона нужно проинтерпретировать значимые зависимости. При наличие в данных пропусков обратите внимание на выбор между casewise and pairwise MD deletion (в чем разница, какие недостатки и достоинства у этих вариантов?).

Учитывая ограниченный объём данных, тем более после разделения по группам (например, наименьшая в private размера 29 без NA), следует воспользоваться pairwise. Более того, из-за признака AVRCOMB использование casewise неявно рассматривает корреляцию в другой подгруппе. Минусом pairwise является то, что из-за попарной проверки могут учитываться данные, пропущенные в других парах.


Затем можно переходить к ранговым коэффициентам корреляции. Расскажите, при каких условиях коэффициенты корреляции Пирсона и Спирмена примерно равны. Приведите примеры, когда один из них больше другого и наоборот. Сравните результаты на ваших данных. Если при сравнении буду найдены заметные различия в результатах, то попробуйте объяснить причину.

compute_corr <- function(group_data) {
  corr_pairwise_pearson <- cor(group_data, use = "pairwise.complete.obs", method = "pearson")
  corr_pairwise_spearman <- cor(group_data, use = "pairwise.complete.obs", method = "spearman")
  return(list(pairwise_pearson = corr_pairwise_pearson, pairwise_spearman = corr_pairwise_spearman))
}

make_corr_plots <- function(corr_list, group) {
  p1 <- ggcorrplot(corr_list$pairwise_pearson,
                   lab = TRUE,
                   lab_size = 1.9,
                   colors = c("blue", "white", "red"),
                   title = paste("Pairwise Cor Pearson Matrix for Group", group),
                   ggtheme = theme_minimal() +
                     theme(plot.title = element_text(size = 6)) # Reduced title size
  )
  p2 <- ggcorrplot(corr_list$pairwise_spearman,
                   lab = TRUE,
                   lab_size = 1.9,
                   colors = c("blue", "white", "red"),
                   title = paste("Pairwise Cor Spearman Matrix for Group", group),
                   ggtheme = theme_minimal() +
                     theme(plot.title = element_text(size = 6)) # Reduced title size
  )
  grid.arrange(p1, p2, ncol = 2, widths = c(4, 4),  # Adjust widths here
               padding = unit(1, "cm"))  # Add padding for extra space
}
library(ggcorrplot)
for (group in names(col_split_mixed)) {
    group_data <- col_split_mixed[[group]]
    cor_matrix <- compute_corr(group_data)
    plot_pearson_spearman <- make_corr_plots(cor_matrix, group)
    print(plot_pearson_spearman)
}

## TableGrob (1 x 2) "arrange": 2 grobs
##   z     cells    name           grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]

## TableGrob (1 x 2) "arrange": 2 grobs
##   z     cells    name           grob
## 1 1 (1-1,1-1) arrange gtable[layout]
## 2 2 (1-1,2-2) arrange gtable[layout]

Теоретическое объяснение (сравнение коэффициента корреляции Пирсона и Спирмена):

В тех ячейках, где совпадают или близки коэффициенты корреляции Пирсона и Спирмена почти нет выбросов и у scatterplot эллиптическая структура, если Пирсон значительно больше Спирмена - есть заметные выбросы. Спирмен может быть больше Пирсона, если выбросы имеются, но их положение не очень значимо относительно других наблюдений.

Private:

  • Выраженный блок GRADUAT, PH_D, SF_RATIO, AVRCOMB, NEW10, OUT_STAT
  • Блок OUT_STAT, R_B_COST. Почему Pearson мощнее Spearman? Наличие нескольких заметных выбросов.

Public:

  • Блок NEW10, AVRCOMB, PH_D, GRADUAT
  • Пара NEW10, R_B_COST

Проинтерпретируйте найденные корреляции - можно ли сказать, что является причиной, что следствием. Если есть какая-то другая причина, которая влияет одновременно на оба признака (скрытый фактор), то попробуйте убрать его влияние с помощью частных корреляций.

Private:

R_B_COST, SF_RATIO (hidden OUT_STAT)

Для R_B_COST и SF_RATIO выберем следующие скрытые факторы:

  • OUT_STAT - стоимость обучения связана с ценой проживания R_B_COST и качеством образования, что влияет на SF_RATIO.

Поставим уровень значимости alpha = 0.05

res<-pcor(col_split_mixed$private[,c("R_B_COST", "SF_RATIO", "OUT_STAT")])
res$estimate["R_B_COST", "SF_RATIO"]
## [1] -0.01433506
res$p.value["R_B_COST", "SF_RATIO"]
## [1] 0.9196585

estimate - получившаяся частная корреляция. Замечаем, что без OUT_STAT коэффициент корреляции Пирсона близок к 0, а p-value ~ 0.90.

AVRCOMB/BOOK (hidden OUT_STAT, NEW10)

  • Для AVRCOMB/BOOK можно выделить OUT_STAT по тем же причинам, что и в предыдущем примере и NEW10, так как они влияют на AVRCOMB
group_1<-c("AVRCOMB", "BOOK")
res<-pcor(na.omit(col_split_mixed$private[,c(group_1, "NEW10")]))
res
## $estimate
##           AVRCOMB        BOOK       NEW10
## AVRCOMB 1.0000000  0.15177454  0.95475775
## BOOK    0.1517745  1.00000000 -0.09376246
## NEW10   0.9547577 -0.09376246  1.00000000
## 
## $p.value
##              AVRCOMB      BOOK        NEW10
## AVRCOMB 0.000000e+00 0.4498169 1.122411e-14
## BOOK    4.498169e-01 0.0000000 6.418045e-01
## NEW10   1.122411e-14 0.6418045 0.000000e+00
## 
## $statistic
##            AVRCOMB       BOOK      NEW10
## AVRCOMB  0.0000000  0.7677672 16.0525714
## BOOK     0.7677672  0.0000000 -0.4708868
## NEW10   16.0525714 -0.4708868  0.0000000
## 
## $n
## [1] 28
## 
## $gp
## [1] 1
## 
## $method
## [1] "pearson"
res_2<-pcor(na.omit(col_split_mixed$private[,c(group_1, "OUT_STAT", "NEW10")]))
res_2
## $estimate
##            AVRCOMB        BOOK   OUT_STAT      NEW10
## AVRCOMB  1.0000000  0.23482510  0.2604398 0.81612411
## BOOK     0.2348251  1.00000000 -0.3815662 0.02045576
## OUT_STAT 0.2604398 -0.38156623  1.0000000 0.27217396
## NEW10    0.8161241  0.02045576  0.2721740 1.00000000
## 
## $p.value
##               AVRCOMB       BOOK   OUT_STAT        NEW10
## AVRCOMB  0.000000e+00 0.24820290 0.19879665 3.728780e-07
## BOOK     2.482029e-01 0.00000000 0.05442782 9.209917e-01
## OUT_STAT 1.987966e-01 0.05442782 0.00000000 1.785854e-01
## NEW10    3.728780e-07 0.92099170 0.17858544 0.000000e+00
## 
## $statistic
##           AVRCOMB       BOOK  OUT_STAT     NEW10
## AVRCOMB  0.000000  1.1834967  1.321494 6.9187347
## BOOK     1.183497  0.0000000 -2.022288 0.1002333
## OUT_STAT 1.321494 -2.0222884  0.000000 1.3856870
## NEW10    6.918735  0.1002333  1.385687 0.0000000
## 
## $n
## [1] 28
## 
## $gp
## [1] 2
## 
## $method
## [1] "pearson"

Получаем, что частная корреляция AVRCOMB/BOOK больше зависит от NEW10, чем от NEW10/OUT_STAT.

Public:

GRADUAT, R_B_COST (hidden OUT_STAT, NEW10)

  • Зависимость между GRADUAT и R_B_COST. Скрытые факторы: OUT_STAT и NEW10.
group_2<-c("GRADUAT", "R_B_COST")
res_1<-pcor(na.omit(col_split_mixed$public[,c(group_2, "OUT_STAT", "NEW10")]), method = "spearman")
res_1
## $estimate
##            GRADUAT  R_B_COST   OUT_STAT      NEW10
## GRADUAT  1.0000000 0.1228687  0.4375959  0.4975827
## R_B_COST 0.1228687 1.0000000  0.2919737  0.2007210
## OUT_STAT 0.4375959 0.2919737  1.0000000 -0.2953932
## NEW10    0.4975827 0.2007210 -0.2953932  1.0000000
## 
## $p.value
##               GRADUAT    R_B_COST     OUT_STAT        NEW10
## GRADUAT  0.000000e+00 0.216289869 3.797030e-06 8.929736e-08
## R_B_COST 2.162899e-01 0.000000000 2.766667e-03 4.205514e-02
## OUT_STAT 3.797030e-06 0.002766667 0.000000e+00 2.451802e-03
## NEW10    8.929736e-08 0.042055136 2.451802e-03 0.000000e+00
## 
## $statistic
##           GRADUAT R_B_COST  OUT_STAT     NEW10
## GRADUAT  0.000000 1.244243  4.890929  5.764986
## R_B_COST 1.244243 0.000000  3.067983  2.059128
## OUT_STAT 4.890929 3.067983  0.000000 -3.107327
## NEW10    5.764986 2.059128 -3.107327  0.000000
## 
## $n
## [1] 105
## 
## $gp
## [1] 2
## 
## $method
## [1] "spearman"

Коэффициент корреляции Спирмена понизился, p-value ~ 0.22. Не отвергаем влияния совокупности скрытых факторов OUT_STAT и NEW10.

OUT_STAT, PH_D (hidden NEW10, SF_RATIO)

  • Коэффициент корреляции Спирмена OUT_STAT и PH_D равен 0.24, хотя визуально кажется, что они коррелируют сильнее. Скрытыми факторами могут быть NEW10, так как NEW10 отражает престижность вузов (NEW10 может зависеть от PH_D, поскольку лучшие студенты хотят найти лучших преподавателей), а OUT_STAT зависит от NEW10 и PH_D.
group_3<-c("PH_D", "OUT_STAT")
res_3<-pcor(na.omit(col_split_mixed$public[,c(group_3, "NEW10")]), method = "spearman")
res_3
## $estimate
##               PH_D    OUT_STAT       NEW10
## PH_D     1.0000000  0.16610803  0.50812933
## OUT_STAT 0.1661080  1.00000000 -0.04604289
## NEW10    0.5081293 -0.04604289  1.00000000
## 
## $p.value
##                  PH_D   OUT_STAT        NEW10
## PH_D     0.000000e+00 0.09520736 4.987348e-08
## OUT_STAT 9.520736e-02 0.00000000 6.458585e-01
## NEW10    4.987348e-08 0.64585848 0.000000e+00
## 
## $statistic
##              PH_D   OUT_STAT      NEW10
## PH_D     0.000000  1.6844818  5.8996943
## OUT_STAT 1.684482  0.0000000 -0.4609178
## NEW10    5.899694 -0.4609178  0.0000000
## 
## $n
## [1] 103
## 
## $gp
## [1] 1
## 
## $method
## [1] "spearman"

Не отвергаем влияние NEW10 на PH_D и OUT_STAT. Можно также рассмотреть частную корреляцию от SF_RATIO, так как это отражает (неявно) процент PH_D и стоимость OUT_STAT за счёт студентов и профессоров.

res_4<-pcor(na.omit(col_split_mixed$public[,c(group_3, "SF_RATIO")]), method = "spearman")
res_4
## $estimate
##                PH_D    OUT_STAT    SF_RATIO
## PH_D      1.0000000  0.20532911 -0.27979699
## OUT_STAT  0.2053291  1.00000000 -0.09968463
## SF_RATIO -0.2797970 -0.09968463  1.00000000
## 
## $p.value
##                 PH_D   OUT_STAT    SF_RATIO
## PH_D     0.000000000 0.02702694 0.002351145
## OUT_STAT 0.027026942 0.00000000 0.287028174
## SF_RATIO 0.002351145 0.28702817 0.000000000
## 
## $statistic
##               PH_D  OUT_STAT  SF_RATIO
## PH_D      0.000000  2.240044 -3.111698
## OUT_STAT  2.240044  0.000000 -1.069669
## SF_RATIO -3.111698 -1.069669  0.000000
## 
## $n
## [1] 117
## 
## $gp
## [1] 1
## 
## $method
## [1] "spearman"

Отвергаем влияние SF_RATIO как скрытого фактора. Не отвергаем, что PH_D скрытый фактор относительно SF_RATIO/OUT_STAT.